Recall from lab that exponentially growing virus, \(V\) should look like this: \[ V(t) = V_0 e^{rt} \] which, on the log scale, is a nice linear model: \[ \log[V(t)] = \log(V_0) + rt \] Our study:
Need to simulate data where \(r\) and \(V_0\) are correlated for each strain.
r_mu <- 1.5 # mean growth rate
r_sig <- 0.5 # sd of growth rate
V0_mu <- 4.6 # mean initial virus population = log(100)
V0_sig <- 1.5 # sd of initial virus population
rho <- -0.65 # correlation between V0 and r
# vector of means
Mu <- c(V0_mu, r_mu)
# vector of standard deviations
sigmas <- c(V0_sig, r_sig) [,1] [,2]
[1,] 1.00 -0.65
[2,] -0.65 1.00
[,1] [,2]
[1,] 2.2500 -0.4875
[2,] -0.4875 0.2500
\[ \begin{bmatrix} a & b \\ c & d \end{bmatrix} \times \begin{bmatrix} e & f \\ g & h \end{bmatrix} = \begin{bmatrix} ae +bg & af + bh \\ ce + dg & cf + dh \end{bmatrix} \]
Simulate parameters by group
nS <- 20 # number of strains
params <- MASS::mvrnorm(n=nS, mu=Mu, Sigma=Sigma)
colnames(params) <- c("V0", "r")
(params <- as.data.frame(params)) V0 r
1 5.6997115 0.8928427
2 4.4961836 1.4243006
3 6.3765896 1.1834381
4 4.9737238 1.9757917
5 7.4924464 1.2853827
6 4.7185409 1.1903997
7 3.4052933 1.9689721
8 2.7186144 1.4827463
9 6.1552155 0.7785731
10 5.1505684 1.3216898
11 6.0818282 1.5490029
12 4.4921862 1.2598660
13 4.3676307 1.3278242
14 0.7001111 2.2344075
15 5.0045218 1.4107406
16 5.9830801 1.4213032
17 3.8837860 1.5354338
18 4.3212154 1.0806057
19 5.1253671 2.1950777
20 3.3916591 1.6799052
Finally time to simulate observations!
Our research questions are:
How would you analyze these data?
a & b vary by clusters, but are independent \[ \begin{align} \log(V_i) &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= a_{\text{Strain}[i]} + b_{\text{Strain}[i]]}\times \text{time} \\ a_{\text{Strain}} &\sim \text{Normal}(\mu_{a},\sigma_a) \\ b_{\text{Strain}} &\sim \text{Normal}(\mu_{b}, \sigma_b) \\ \mu_{a} &\sim \text{Normal}(4,1.5) \\ \sigma_{a} &\sim \text{Exponential}(2) \\ \mu_{b} &\sim \text{Normal}(1,1) \\ \sigma_{b} &\sim \text{Exponential}(3) \\ \sigma &\sim \text{Exponential}(2) \end{align} \]
a & b vary by clusters, but are correlated $$ \[\begin{align} \log(V_i) &\sim \text{Normal}(\mu_i, \sigma) \\ \mu_i &= a_{\text{Strain}[i]} + b_{\text{Strain}[i]]}\times \text{time} \\ \left[\begin{matrix} a_{\text{Strain}} \\ b_{\text{Strain}} \end{matrix} \right] &\sim \text{MVNormal}\left(\left[\begin{matrix} \mu_{a} \\ \mu_{b} \end{matrix} \right],\Sigma \right) \\ \Sigma &= \left( \begin{matrix} \sigma_a & 0 \\ 0 & \sigma_b \end{matrix} \right) \text{Rho} \left( \begin{matrix} \sigma_a & 0 \\ 0 & \sigma_b \end{matrix} \right) \\ \mu_{a} &\sim \text{Normal}(4,1.5) \\ \sigma_{a} &\sim \text{Exponential}(2) \\ \mu_{b} &\sim \text{Normal}(1,1) \\ \sigma_{b} &\sim \text{Exponential}(3) \\ \text{Rho} &\sim \text{LKJcorr}(2) \\ \sigma &\sim \text{Exponential}(2) \end{align}\] $$
ulam()a & b vary by clusters, but are independent
ulam()a & b vary by clusters, but are independent
mean sd 5.5% 94.5% n_eff Rhat4
a[1] 5.6197359 0.28164215 5.1356684 6.0741868 449.9927 0.9981760
a[2] 4.7118293 0.28615075 4.2662249 5.1606819 487.9817 0.9980835
a[3] 6.2525709 0.26581581 5.8259082 6.6684760 375.2530 0.9980828
a[4] 4.6659927 0.27625705 4.2391181 5.1218754 445.6447 0.9982450
a[5] 7.4975931 0.28200501 7.0029107 7.9305953 446.4993 1.0010159
a[6] 4.6477172 0.28387578 4.2254237 5.1334320 515.7942 0.9989191
a[7] 3.8892686 0.28299326 3.4693857 4.3807230 557.8814 1.0045932
a[8] 2.9734650 0.27109802 2.5508475 3.3999261 376.7245 1.0004947
a[9] 5.6376343 0.29488779 5.1884568 6.0771947 446.8889 1.0029273
a[10] 5.1894291 0.27536434 4.7141119 5.6159996 438.3649 1.0018388
a[11] 5.8634807 0.27772760 5.3923784 6.3354943 451.7656 0.9982016
a[12] 4.1641750 0.29046496 3.7086837 4.6194622 470.7299 0.9980760
a[13] 4.3872350 0.27518780 3.9770622 4.8469286 505.0801 0.9989898
a[14] 1.2771100 0.30963280 0.8441708 1.7853930 307.5220 0.9998305
a[15] 5.2451792 0.25533451 4.8367318 5.6186080 401.0076 0.9983295
a[16] 5.8120824 0.26097571 5.3944576 6.2242712 642.1812 0.9980392
a[17] 4.0266713 0.28655454 3.5669885 4.5023922 467.8292 0.9994088
a[18] 4.2720804 0.26584935 3.8388472 4.6814475 404.6871 0.9988240
a[19] 5.7211973 0.27967798 5.2838661 6.1614799 610.1397 0.9992123
a[20] 3.5095258 0.28534628 3.0684600 3.9735664 372.2920 1.0034292
a_mu 4.7252093 0.29493601 4.2705324 5.2030886 697.2526 0.9988166
a_sd 1.3396701 0.20974641 1.0563460 1.6870729 508.1255 1.0015487
b[1] 1.0294046 0.08191225 0.8973560 1.1581166 409.2617 0.9993012
b[2] 1.4568568 0.08550066 1.3160524 1.5857219 499.0913 0.9980833
b[3] 1.2596826 0.08064037 1.1295799 1.3928521 370.2664 0.9984688
b[4] 2.0746709 0.07813137 1.9473135 2.1932221 372.4381 0.9987550
b[5] 1.2281400 0.08432477 1.0976268 1.3665780 364.2279 1.0003461
b[6] 1.2463987 0.08155160 1.1141551 1.3718674 494.7723 0.9996309
b[7] 1.8302930 0.08344857 1.6939483 1.9615252 504.1925 1.0054988
b[8] 1.3686721 0.07815777 1.2420503 1.4925317 372.4234 1.0021726
b[9] 0.8990466 0.08866126 0.7627035 1.0444876 526.6238 1.0018244
b[10] 1.3045427 0.08006247 1.1831102 1.4277186 468.2089 0.9998287
b[11] 1.5981687 0.08279476 1.4689424 1.7377499 498.9307 0.9983048
b[12] 1.2996811 0.08305029 1.1676456 1.4362483 467.3709 0.9981242
b[13] 1.3551176 0.08227767 1.2229513 1.4774558 505.2874 0.9988714
b[14] 2.0515093 0.09284407 1.9039347 2.1831765 302.6048 1.0009567
b[15] 1.3380622 0.07448518 1.2197519 1.4601463 425.0079 0.9984582
b[16] 1.4086105 0.07839334 1.2891097 1.5310297 680.2872 0.9981991
b[17] 1.4665637 0.08284514 1.3389944 1.5871137 470.3861 0.9988451
b[18] 1.1830804 0.07722528 1.0549612 1.3046221 397.6879 0.9979990
b[19] 1.9731130 0.07944392 1.8428377 2.0924031 509.8644 1.0024707
b[20] 1.6173833 0.08085352 1.4745344 1.7402543 361.8035 1.0021741
b_mu 1.4483017 0.08322592 1.3281445 1.5826809 541.2091 1.0039075
b_sd 0.3454524 0.06676705 0.2545205 0.4594744 279.2551 0.9981472
sigma 0.2301938 0.04110536 0.1732079 0.3027420 114.1453 0.9980733
a & b vary by clusters, but are correlated
mean sd 5.5% 94.5% n_eff Rhat4
b[1] 1.0143803 0.08744104 0.8837770 1.15801625 721.4040 0.9987863
b[2] 1.4569767 0.08179049 1.3351552 1.59187810 536.6146 1.0006891
b[3] 1.2426888 0.08137205 1.1261775 1.37907975 620.6789 0.9994441
b[4] 2.0800591 0.08071669 1.9575744 2.22099455 650.7942 0.9992440
b[5] 1.2136164 0.07634753 1.0985469 1.32535325 629.5236 0.9999293
b[6] 1.2392340 0.07849195 1.1139775 1.35824845 486.3734 0.9988040
b[7] 1.8445091 0.08234411 1.7169503 1.97709040 704.2020 1.0004837
b[8] 1.3839925 0.07384597 1.2634173 1.49657395 577.2652 0.9980004
b[9] 0.8851179 0.07599035 0.7585436 1.00713860 487.2414 0.9980013
b[10] 1.3015516 0.07936911 1.1725432 1.43605375 598.3270 0.9985473
b[11] 1.5967882 0.07801285 1.4753589 1.71759745 527.3257 1.0000229
b[12] 1.3017448 0.07661439 1.1846277 1.42124360 673.8910 0.9984481
b[13] 1.3535050 0.07775585 1.2292371 1.47398155 657.0217 0.9988271
b[14] 2.0786186 0.08421506 1.9410978 2.20668740 472.8595 0.9987464
b[15] 1.3379790 0.07069365 1.2325435 1.45309235 529.5274 0.9980077
b[16] 1.4011446 0.07276654 1.2914682 1.51331075 649.7329 0.9980513
b[17] 1.4734786 0.07084068 1.3650928 1.58876935 750.1340 0.9985369
b[18] 1.1817457 0.07258477 1.0719075 1.29910425 484.5213 0.9980433
b[19] 1.9719898 0.08109897 1.8410768 2.09867885 1118.9537 0.9992206
b[20] 1.6216965 0.07526039 1.5038094 1.73946045 543.1927 0.9990044
a[1] 5.6774484 0.29440941 5.2285598 6.11345725 766.1885 0.9981898
a[2] 4.6985481 0.27857557 4.2427601 5.13432750 601.9208 1.0005713
a[3] 6.2984875 0.26145359 5.8670361 6.69775910 613.5075 0.9989312
a[4] 4.6387087 0.26812261 4.1674365 5.06811355 550.1051 0.9990865
a[5] 7.5493572 0.26083686 7.1196535 7.94289415 692.7991 0.9984522
a[6] 4.6603255 0.26946141 4.2316747 5.09041660 518.9945 0.9996686
a[7] 3.8441048 0.28784602 3.3698946 4.29179275 759.3335 0.9992679
a[8] 2.9212086 0.25510791 2.5156858 3.31248930 577.8133 0.9980125
a[9] 5.6906856 0.26501390 5.2771532 6.10849150 520.9542 0.9981407
a[10] 5.1970856 0.28115695 4.7533679 5.63334530 671.4429 0.9981073
a[11] 5.8557941 0.26654734 5.4458153 6.27611195 551.3730 1.0046470
a[12] 4.1716393 0.26398835 3.7553513 4.57763940 678.5736 0.9980023
a[13] 4.4036094 0.25606265 3.9938437 4.83258585 623.0677 0.9980008
a[14] 1.1732646 0.28256899 0.7603441 1.63746960 459.9166 0.9979982
a[15] 5.2541280 0.24457845 4.8640175 5.66460605 516.2838 1.0005499
a[16] 5.8344220 0.25412395 5.4427746 6.24292575 749.5599 0.9990374
a[17] 4.0016485 0.24891853 3.6205492 4.38283615 802.1613 0.9985174
a[18] 4.2647436 0.26098545 3.8904597 4.68331845 624.8952 0.9985104
a[19] 5.7149885 0.29141917 5.2949689 6.20065335 1091.2689 0.9982974
a[20] 3.5101416 0.26644409 3.0791245 3.95307420 584.6252 0.9983197
a_mu 4.7318570 0.30892371 4.2267149 5.20004140 777.5271 0.9981810
b_mu 1.4486811 0.07625235 1.3275034 1.57009370 640.7886 0.9986301
sigma_ID[1] 1.3803440 0.21540617 1.0699386 1.73998950 389.9947 0.9980070
sigma_ID[2] 0.3493480 0.05968775 0.2621147 0.45054538 483.9477 0.9986399
Rho[1,1] 1.0000000 0.00000000 1.0000000 1.00000000 NaN NaN
Rho[1,2] -0.3844072 0.17204899 -0.6324356 -0.07074392 491.9179 0.9988999
Rho[2,1] -0.3844072 0.17204899 -0.6324356 -0.07074392 491.9179 0.9988999
Rho[2,2] 1.0000000 0.00000000 1.0000000 1.00000000 NaN NaN
sigma 0.2233401 0.03795678 0.1694872 0.29153988 154.1904 1.0168633
red is estimate without correlation blue is estimate with correlation
m3 <- ulam(
alist(
V ~ dnorm(mu, sigma),
mu <- (a_mu + alpha[ID, 1]) + (b_mu + alpha[ID, 2])*time,
# priors
# adaptive priors - non-centered
transpars> matrix[ID, 2]:alpha <-
compose_noncentered( sigma_ID, L_Rho_ID, z_ID),
matrix[2,ID]:z_ID ~ normal( 0 , 1 ),
a_mu ~ dnorm(4,1.5),
b_mu ~ dnorm(1,1),
vector[2]:sigma_ID ~ dexp(2),
cholesky_factor_corr[2]:L_Rho_ID ~ lkj_corr_cholesky(2),
sigma ~ dexp(2),
# compute ordinary correlation matrices from Cholesky factors
gq> matrix[2,2]:Rho_ID <<- Chol_to_Corr(L_Rho_ID)
), data = df
) mean sd 5.5% 94.5% n_eff Rhat4
Rho[1,1] 1.0000000 0.000000 1.0000000 1.00000000 NaN NaN
Rho[1,2] -0.3844072 0.172049 -0.6324356 -0.07074392 491.9179 0.9988999
Rho[2,1] -0.3844072 0.172049 -0.6324356 -0.07074392 491.9179 0.9988999
Rho[2,2] 1.0000000 0.000000 1.0000000 1.00000000 NaN NaN
mean sd 5.5% 94.5% n_eff Rhat4
Rho_ID[1,1] 1.0000000 0.0000000 1.0000000 1.0000000 NaN NaN
Rho_ID[1,2] -0.3795035 0.2037607 -0.6572971 0.0284156 122.0967 1.042638
Rho_ID[2,1] -0.3795035 0.2037607 -0.6572971 0.0284156 122.0967 1.042638
Rho_ID[2,2] 1.0000000 0.0000000 1.0000000 1.0000000 NaN NaN